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SUMMARY 


An analytical investigation of the vibrations of thermally- 
stressed orthotropic plates in the prebuckled region is 
presented. The investigation covers the broad class of trape- 
zoidal plates with two opposite sides parallel. Each edge 
of the plate may be subjected to different uniform boundary 
conditions. Variable thickness and arbitrary temperature 
distributions (analytical or experimental) may be prescribed. 
Generality is achieved in the analysis through the treatment 
of boundary conditions, the choice of functions for stress 
distributions and deflection distributions, and the use of 
numerical integration for the evaluation of matrix elements. 
Results obtained using this analysis are compared to experi- 
mental results obtained for isotropic plates with thermal 
stress, and to results contained in the literature for 
orthotropic plates without thermal stress. Good agreement 
exists for both sets of comparisons. Calculations for 
several orthotropic plates with thermal stresses indicates 
that the effect of orthotropy on the frequencies may be 
large and should not be ignored. 
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LIST OP SYMBOLS 


- stress matrix from complementary 
energy 

- elements of [A] 

- aspect ratio, length squared/area 

- plate length 

- matrix of material constants in equa- 
tion for strains in terms of stresses 

- elements of La], Equation (3) 

- generalized stiffness matrix from 
bending energy 

- element of [B] 

- plate dimension measured along left 
edge from x-axis to. top corner 

- plate dimension, plate width at left 
edge minus b 

- ratio of b^/b 

1] h 

- b^ (x,y) = 0, equation of i por- 
tion of plate boundary 

jr 

- coefficient of the pq term of the 
assumed stress function solution 

- matrix of material constants in 
equation for stresses in terms of 
strains, [fi] =• [a]" 1 

- elements of [E] 

- stress function solution to the inplane 
e q u ilibri urn e qua 1 1 o n s 

- frequency of i ti:i mode, cycles per 
second 


yii 



f 


g 



h 


h 


r 


CM] 


M. . , , 

XJ ,k l 


- f(x,y) - function which forces the 
assumed stress function solution to 
satisfy the stress boundary condi- 
tions 

- g(x,y ) - forces the assumed dis- 
placement solution to satisfy the 
displacement boundary conditions 

th 

- coefficient of the ij term in the 
assumed displacement function solution 

- h(x,y) - function to represent any 
variation in plate thickness 

- plate thickness at some reference 
point 

- mid-plane energy matrix, associated 
with the thermal stresses moving 
through small out-of-plane displace- 
ments 

- elements of [M] 

h h h 

2 2 2 

- / ^ CT-^dz ,/• ^ o^dz,f T i2 dz res P ectl vely 

~2 ~2 ~2 


n 


[T] 

T 

ij ,k l 
T 


1 re f 


cr 


- coordinate normal to plate boundary 
(in-plane) 

- generalized mass matrix from kinetic 
energy 

- element of [T] 

- T(x,y) - difference between the tem- 
perature at a point (x,y) on the 
plate and the original, uniform ref- 
erence temperature (non-dimensional) 

- difference between the temperature 
at some reference point on the plate 
and the original, uniform reference 
temperature 

- The magnitude of T pe f at which the 
free vibration frequency vanishes. 

By definition, the. thermal buckling 
temperature. 
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time 


U' , V , W 
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ij 
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AT 


£ 1 5 S 2 
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- displacements in the x,y, and z 
directions respectively 

- independent in-plane variables 

- angle between plate leading edge 
and x-axis , measured positive 
counter- cl o c kwis e 

- taper parameter, a = ^ tan a 

- coefficient of thermal expansion in 
x arid y directions respectively 

t/ 

- i j term of general assumed dis- 

placement function 

- angle between x-axis and tire line 
dividing the plate for thickness 
distribution purposes 

- non-dimensional form of 3 

- thermal loading matrix 

- element of {!} 

- angle between plate trailing edge 
and x-axis , measured positive 

c o unte r-clockwise 

- taper parameter, y - tan y 

- shear strain 

- pq term of general assumed stress 
function 

- an increment of T _ , gives magni- 
tude of the temperature distribution 
under consideration 

- normal strains in the x and y direc- 
tions respectively 

- non-dimensional independent space 
variable, n = y/b 
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- non-dimensional independent space 
variable, £ = x/a 

- energy of a system per unit t imo 

- complimentary energy 

- plate material density, mass per 
unit volume 

- normal stresses in x and y directions 
respectively 

•- shear stress 


-- vibration frequency of the thermally 
stressed plate, radians/sec. 

vibration frequency of the plate at 
T = 0 , radians/sec . 


MOTE ON SUBSCRIPT CONVENTION 

Numeric subscripts indicate the component of a quantity 
in a coordinate direction (e . g. , - normal stress in the 

1 or x - direction). A subscript of x,y,£, or n denotes 
dif ferentiation.gwith respect to that independent variable 
(e.g. , (a^) = ( a i))* All other alphabetic subscripts 

(i, j, k, , p, q. etc.) will refer to either terms in a 
series or elements in a matrix. 
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I. INTRODUCTION 


Considerable work has been reported in the literature 
on the problem of finding the frequencies and modes of vi-" 
bration of a rectangular orthotropic plate at ambient 
temperature. A combination of the work of Hearmon (Ref. 1, 

2, 3); Hoppman, Huffington, and Magness (in various combina- 
tions, Ref. 4, 5, 6, 7 , 8); Kanazawa and Kawai (Ref. 9) and 
Wah (Ref. 10) provide solutions for rectangular plates with 
any boundary condition except completely free. 

In contrast, no literature was found pertaining to the 
free vibration frequencies of an orthotropic plate subjected 
to a thermal loading. For the special case of a thermally 
stressed isotropic plate, the torsion mode of the plate with 
cantilever boundary conditions has been rather thoroughly 
investigated (Ref. 11, 12, 13, 14, 15). 

Ref. 16 presents an analysis of thermally stressed iso- 
tropic plates for various boundary conditons , ranging from 
plates completely clamped through several combinations of 
mixed boundary conditions to plates with all edges completely 
free. This paper extends the analysis of Ref. 1 6 to include 
orthotropic plates with a thorough discussion of the associated 
computer program. 


In the sense that both computability and equilibrium 
are satisfied as closely as one pleases at every point inter- 
ior to the plate and on the boundary , this paper presents an 
analysis that provides, in a practical computational sens.e, 
a solution to the thermally stressed plate vibration problem 
for all trapezoidal plates with two opposite sides parallel, 
and with one of the axes of elastic symmetry parallel to, 
these sides, restrictions that could be. easily relaxed. 

The analysis and associated computer program are of 
sufficient generality that isotropic plates are included as 
a special case of orthotropic plates. Various boundary 
conditions may be arbitrarily assigned to the different sides 
of the plate. Thus, the solution for the vibrations of 
thermally stressed plates with boundary conditions ranging 
from completely clamped to completely free with any 
combination of clamped, pinned, and/or free edges may be 
obtained. 

A small number of quantitative strain measurements (not 
included herein) plus the abundance of experimental dynamic 
response data for various pi an form shapes and boundary condi- 
tions of isotropic plates indicates that the' stress distributions 


1 



as determined herein are correct. 


No thermally stressed orthotropic plate data, either 
analytical or experimental, were found in the literature; 
however, for orthotropic plates without thermal stress, 
comparison is made to both analytical and experimental data 
from the literature. Further comparison is made with experi- 
mental data for several modes of thermally stressed isotropic 
plates with various planform shapes, boundary conditions, and 
temperature distributions. 


II. APPLICATION OF THE ENERGY EQUATION 


Because of a difference in notation between that used 
herein and that used in other sources, the derivation of the 
expressions for potential and complementary energy will be 
shown. Except for this difference, the procedures used are 
well known. Additional details, if desired, may be found in 
Ref. 17- 

A. Potential Energy 

The forces are taken as the independent variables and 
the variation of the total energy is taken with respect to 
the displacements. With the assumptions of plane stress 
and no body or surface forces. 


6 TT 


///{(o 1 Se 1 + a 2 Se 2 


+ t 12 «y 12 ) 


+ p 


w <$w} dxdydz (1) 


The orthotropic stress-strain relations will be taken 

as , 

{ £ } = [a] {a} + {a} T (2) 


where 




The inverse of eq. (2) is 

{a} = [E] {e} - [E] { e}T (4) 

The Von Karman strain-displacement equations are used. 



Substitute eqs . 4 and 5 into eq. 1, carry out the indicated 
operations, neglect fourth order terms, and integrate through 
the thickness to get. 
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Sit = yr o//{N u + N v + N (u + v ) 
2 ix 2 y 1 2 y x 


+ ~ [E (w ) + E (w ) + 2E w w + 4G (w ) 

12 ii xx 22 yy 12 XX yy iz K xy' 
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+ 2 ph w w} dxdy = 0 
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N = 
i 


h 

2 

h 

'2 


[E u + E v - T (E a -+ E a ) ] dz 

11X 12 y 111 122 
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h 

2 
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[E u + E v 
12 X 22 y 


T (E a + E a )] dz 
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(7) 


N = 

1 2 
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2 

h 

'2 


G i 2 (u y + v x> dz > 


Taking the variation with respect to u gives. 


ff ( H Su + N <5u ) dxdy = 0. 
l x 12 y 


Integrating this result by parts, noting that for any 
solution to a particular problem the boundary conditions must 
be satisfied, leaves the in-plane equilibrium equation in the 
x-direction. 



(8b) 


Performing a similar series of operations on v gives 
for the y-direction. 


(N ) + 

1 2 X 


(N ) 

2 y 
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These eqs. (8) have the solution, F, such that. 


N 

i 




( 9 ) 


N = 
1 2 




Thus the variational expression for potential energy of 
a thermally stressed, orthotropic plate becomes. 


= | SJflgr [E u (« n ) ! + E 2! (w yy ) 2 + 2E i2 W xx w yy 


+ 1 ,a i2 (w xy ) 3 


( 10 ) 


2 2 

+ [ P (w ) + P (w ) — 2F w w 1 

L yy x x' xx y xy x y J 


+ 2 ph w w } dxdy = 0 


B. Complementary Energy 

In developing the expression for complementary energy, the 
unknown forces are varied and the displacements are held 
constant. Thus, with the same assumptions as used in the 
treatment of potential energy, 

<5tt = ///{.e So + e So + y 6x Idxdydz - 6w n = 

1 1 2 2 12 12 B 


0 



where Wg represents the work done by the stresses on the 
portion of the boundary on which the displacements are speci- 
fied. In this treatment, if the displacements on any part of 
the boundary of the plate are to be specified, they will be 
specified to be zero. Thus Wg = 0 . 

Substitute eq. ( 4 ) for the stresses to get. 


Sir* = fffie 6 [E e + E e - (E ct + E a )T] 
1111 122 111 122 

+ e <5 [E £ + E £ - (E a + E a )T] 

2 222 121 121 222 

+ y 5 G y } dxdydz - 0 
12 12 12 


Because the bending stresses have been expressed in terms 
of the displacement only and small deflections are assumed, 
the stresses that remain in the equations are not functions 
of the out-of-plane displacements; they are "membrane stresses" 
resulting only from the in-plane displacements and/or tempera- 
ture. Thus, only the linear strains resulting from in-plane 
deformation need to be considered and eq. ( 5 ) can be simplified 
to , 


e i " u x 


e = v 

2 y 


= (u y + V 


With these relations, the complementary energy can be 
expressed as. 


< 5 it = ///{u 6 [E u + E v - (E a + E a )T] 

X 1 1 X 12 y 111 122 


+ V 6[E u v + E v„ - (E a + E a )T] 
y 12 X 22 y 121 222 


+ (u y + V x ) «[a i2 (u y + v x n 


dxdydz - 0 



Integrate through the thickness, substitute equations 
(7) and. (9) and define the strains to be. 


e = u = ;~(a P +a P ) + a T 
1 x h i i yy 12 xx 7 i 


e = v t = r- (a F „ + a P )+aT 
2 y h 22 xx i2 yy 2 


Y i2 = u y + v x = 


h b i Z F xy 5 


from which the complementary energy for a thermally stressed, 
orthotropic plate becomes. 


Sir* = 677 { ir- [a (F)* + a (F)* + 2 a F F vv 
2 h 11 yy 22 xx 12 xx yy 


+ b ,a (F xy) ^ 


( 11 ) 


+ (a F + a P ) T } dxdy = 0 
1 yy 2 xx 


C. The Equations in Matrix Form 

Consider first the potential energy, eq. (10). Assume a 
displacement function of the form, 

N M 

w(x,y,t) - E E h. (t) a.. (x,y) (12) 

i=o j-o * J 

where each a.. (x,y), (1) satisfies the displacement boundary 
conditions, pi?) is continuous, and (3) has at least continuous 
first derivatives. 

Substitute this into eq. (10), take the variation with 
respect to h k £ and collect coefficients of like h^ • to get the 
matrix equation. 
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[B] {h ± .} + K^M] {h ij .} - X 2 [T] ih ± .} = 0 (13) 


where the non-dimensionalized matrix elements and associated 
parameters are given in Appendix A. 

Now assume a stress function of the form, 

F (X > y) = j„ J„ C pq T pq < X ’ y) <X1|) 

p=o q=o ^ * 1 


where each y Da (x,y), (1) satisfies the stress boundary condi- 
tions, (2) is q cont inuous , and (3) has at least continuous 
first derivatives. 

Substitute this into eq. (11), take the variation with 
respect to C rs and collect coefficients of like C p q to get 
the matrix equation. 


[A] {C pq } + K 2 m = 0 


(15) 


where the matrix elements are also given in Appendix A. 

Thus, given a temperature distribution, {T} can be 
calculated, eq. (15) can be solved for {C pq }, and values of 
the derivatives of the stress function can 4 be found. Using 
this information, the elements, My can be calculated and 

eq. (13) can be solved for the vibration frequencies and modes 
with thie buckling mode and AT cr obtainable as a limiting case 
when x| - 0 . 

D. Deflection Function and Stress Function 


At this point, a choice will be made concerning the 
form of the assumed deflection function and stress function. 

By observing the physical system, it can be seen that the 
deflected surface of the plate and the stresses within the 
plate are continuous and have at least continuous first 
derivatives. Thus, the functions to be assumed as solutions 
to the problem must belong to the class of functions which 
are continuous and have at least continuous first derivatives. 
The assumed solution must also satisfy the boundary conditions 
discussed in the next section. 
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A truncated power series in the independent space varia- 
bles satisfies the continuity requirements. Thus, the functions 
assumed for the deflection, w, and for the stress function, 

F, will be truncated power series. 

E. Boundary Conditions 

The polynomial resulting from a truncated power series 
will not in general satisfy the boundary conditions. There- 
fore, the polynomial representation must be modified by an 
additional function which forces satisfaction of the required 
boundary conditions. 

Let the displacement function have the form, 

N M . 

w (x,y) = g (x,y) E E h. .x y J 

i=o j=o J 

where g(x,y) is the boundary condition function which insures 
satisfaction of the displacement conditions at the boundary. 

The stress function will have the form, 

S T 

P (x,y) = f (x,y) E E C x p y q 

p=o q=o 


where f(x,y) is the boundary condition function which insures 
satisfaction of equilibrium in the plane of the plate at the 
boundary, i.e., the stress boundary condition. The specific 
form of each will now be considered. 

1. Displacement Function 

Three types of displacement boundary conditions are 

considered herein: 

(a) Both displacement and slope normal to the edge 
of the plate are assumed to be zero; i.e., the 
edge is clamped. 

(b) Only the displacement is assumed to be zero 
and the slope is left unspecified resulting in 
a pinned (simply-supported) condition. 

(c) Both slope and displacement are left unspecified 
leaving the edge completely free. 
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Now, given a particular plate geometry, the equation of 
the boundary may be expressed as a polynomial, say. 


b (x,y) = 0 . 


Therefore, in order to force the displacement to be zero on, 
the boundary, simply let. 


g ( x ,y ) = b (x,y) , 

so that for any point on the boundary, (x R ,y R ), the deflec- 
tion will be a u 

N M . . 

w (x B ,y B ) = g (x B ,y B ) Z ,Z h^ x g y" = 0 


This satisfies condition (b) because the first deriva- 
tive, 1^- , will not in general be zero but will be left to 
take on whatever value is required for a minimum energy con- 
figuration. 

Condition (a) may be satisfied by letting 

2 

g (x,y) = [b(x,y)3 


The displacement will again be zero, but now the first de- 
rivative will also be zero on the boundary: 


I ■ j 0 j 0 *ij + x y ssi&ii 

N M 2 a ^ . 

- Z Z h.. { [ b ( x , y ) ] — (x y J ) + 2 x^bCx^) x 


i-o j-o 


ij 


3ri b(x ^ )} 

and at a point (x g ,y g ) on the plate boundary, 


3w 

9n 


0 
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Case (c) may be satisfied simply by letting 

g(x,y) 21= [b(x,y)]° 


Thus, in the case of an edge free to displace out of the plane 
of the plate, both the displacement and slope will be left to 
take on whatever values are required for a minimum energy 
configuration. 

If the plate is a polygon of N sides, write 


g(x,y) 


N k 

tt [b. (x,y) ] 
i=l 1 


t h 

where b.(x,y) is the equation of the i side of the poly- 
gon. (The sides need not be straight..) kj[ will be either 
0, 1, or 2 as described above. 

2. Stress Function 


Two types of stress boundary conditions are consid- 
ered herein: 

(a) The in-plane stresses normal to a boundary are 
specified to be zero. That is, the plate is 
left free to expand in the in-plane direction. 

(b) The stresses are completely unspecified or, 
equivalently, the in-plane displacements are 
specified to be zero. The stresses will take 

on whatever values are required for satisfaction 
of equilibrium. 

Thus, condition (a) will be termed "free" and (b) will 
be termed "clamped". These conditions are fulfilled in a 
fashion similar to that used with the deflection function. 

Recall from classical elasticity theory that the 
stresses normal to the edge of a plate will be zero if the 
stress function and its first derivative (normal to the 
edge) vanish there. Therefore, on any portion of the boundary 
on which a free condition is desired, the equilibrating 
function is, 

f (x,y ) - [b 1 (x,y) ] = 0 
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where, as before, bj_(x,y) = 0 is the equation of that portion 
of the boundary. 

A clamped condition on any portion of the boundary can 
be satisfied by setting. 


f(x,y) = [b^(x,y ) ]° = 1 . 


Thus, as was done with the boundary condition function, 
the equilibrating function is 

N k. 

f(x,y) = n [b. (x,y) ] 1 
1=1 1 


where kj; = 0, 2 for clamped or free conditions respectively 
on the "ith" side of the polygon. 

With these conditions, it is now possible to specify 
six different types of boundary conditions on any plate edge. 
The first letter of the notation used herein will denote the 
displacement boundary condition by using P = free, P = pinned 
or simply supported, and C = clamped. The second letter 
denotes the stress condition so that designations possible for 
any given edge are : 


DESIGNATION 

DISPLACEMENT 

STRESS 

(1) 

P - P 

free 

free 

(2) 

P - P 

pinned 

free 

(3) 

1 

o 

clamped 

free 

(4) 

p - c 

free 

clamped 

(5) 

p - c 

pinned 

clamped 

(6) 

c - c 

clamped 

clamped 


I.II . PROGRAMMING 


A. The Equations and the Plate Geometry 

The equations to be programed are eqs . 13 and 15 with the 
matrix elements as given in Appendix A. 
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The planform and descriptive parameters of the plates 
considered herein are shown in Pig. 1, The restriction that 
two of the sides are parallel' was made to simplify the numerical 
integration scheme. 

The plate edges are numbered clockwise (in the top view) 
beginning with the edge containing the origin. Thus, the.' 
equations of the four edges as used in the boundary condition 
functions are:- 


b ^ (S,n) = 5 = 0 

b (£,n) = (i + y5 - n) = o 

2 

b (5,n) = (l - 5) = 0 

3 

b (5,n) = (b - y 5 + n) = o 

4 1 


( 16 ) 


B. Logic Plow Diagram 

The organization of the parts of the program is presented 
in a logic flow diagram shown in Appendix B. It should be 
noted first that if several sets of material properties and/or 
aspect ratios are to be investigated, the [B] and [A] matrices 
need not be integrated each time if the integrant are treated 
as four separate terms. Each of these terms need, only be 
integrated once, then multiplied by the appropriate constant 
and added together to make up the whole integral for either 
isotropic or orthotropic materials. The [T] and [M] matrices 
are independant of both the material properties and aspect 
ratio. The program is structured to Include this feature. 

A complete listing of the program is given in Appendix C. 

C. Programming Boundary Conditions 

The subroutine "FUNCTN" which calculates the displacement 
boundary condition function and the equilibrating function and 
their derivatives is very straightforward . Since both func- 
tions have the same form, the same subroutine can be used to 
simplify the user's task of calculating the exponents required. 
There is no need, for example, to remember that a zero exponent 
on the stress function means a clamped edge while the same 
exponent in the displacement boundary condition function means 
a free edge . 
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The first step was to write down the function and its 
five derivatives, leaving the exponents as variables. For 
example, 

d It I. 

F = 5 (1+aC-n) (1-0 (b^-yC+n) , 

3F I ■ - 1 1 3 1 j, 

3 ^- = -i 2 C5 (l+.a5-n) 1 d-O (b^-y^+n) ] 


ii I, I I -1 

+ ICC ( 1 +aC-n) (l-C) 3 (b -yO-n) J 

4 1 


I I -2 I i k 

3 2 F/an 2 = 1 (I - 1 ) U 1 ( l+a£~n ) 2 (1-5) 4 (b -yOn) ] 

2 2 1 


I 1 J 2 I 3 I ,"2 

+ 1(1 -1) [5 1 (l+a5-n') (1-5) (b -y?+n ) 3 

4 4 1 


1 1 1 2~ 1 1 3 1 - 1 

- 21 I [5 (l+a5-n) ( 1 - 5 ) ■ (b -y5+n) 4 3 

2 4 1 


Thus, it can be seen that Ij_, I±-l, and Ij_-2 are re- 
quired for calculating the function and its derivatives. In 
the subroutine, the variable IEX( I , J) contains these quantities. 
The "I" refers to the four factors making up the function and 
"J" to I 3 -O, Ij_-1, or Ii-2. This is done in "DO-LOOP” 
number five. 

Next, this information is used to calculate all the 
factors T (M,K) required for the function and its derivatives. 
For example, 

T (1,1) = 5 1 

I -1 

T (1,2) = 5 1 

I 3 

T (3,1) = (1-0 

1-2 

T (4,3) = (b +y5-n) 
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Finally s this information is used to calculate the function 
and its derivatives. 

D. Numerical Integration 

Integration of the elements of the various matrices is 
performed using the Gaussian Quadrature rule.. The plate is 
divided into two parts by the line at angle 0. This provides 
for more accurate results when the leading or forward part 
of the plate has a different thickness function than does 
the rearward part. Since the Gaussian Quadrature is defined 
on the interval [-I, 1] it is necessary to transform the 
points and coordinates. 


If u 


= <Kx) then f l f (u) du = /V x f O(x) ] dx 


Then if 


1 N 

/_ 1 f ( x ) dx * Z w v f (xj s 


k=l 


k * v "k' 


where 


b 


N 


U f (u)du “^i \ F( V 


\ ‘ Hr 1 V u k - * (x k> 


For this problem. 


lb a 

/ (/ f(?,n)dn + I f (?,n )dn}dc 


11 Zip 

= f {/ f[<S)(x,y) , (x,y) ] -r- 1 - dy 

-1 -1 1 3y 


1 3 ^ . 

+ ./ (x,y) , 'l> 2 (x,y)l dy gf dx , 
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where 


Thus 


a = 35> b = l+a£, c =y£ -b 


5 = 1/2 (x+1) 


( 17 ) 


Hence, 


n = 1/2 {y [ 1 + (a-3)5] + 1 + (a+3)5h ( 35 £h£l+a 5 ) 


n = 1/2 {y[b + (3-y)?]-b + (3+y)) > (y?ln£35) 


9 (j> 

9x 


= 1/2 




x 


3y 


1/2 [l+(a-3)5] 


9 ^ 


3 * • 1/2 Cb^ + (B-y)?] 


and the Integrals may be evaluated by 


f( 5 ,n) dn + / c f(5,ri) d n } d£ 


N w, N w T 

K r r. Lt 


l -f- { Z -f [b + ( 3 -y) 5 k ] fU k ,n L ) 

k=l 2 L=1 2 1 K K L 


2N W M 

+ S -5 [1 + (a-3)5 k l f(5 k ,n M )> 


M=N+1 
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where w k , w L , and wjyj are the values on [-1,1] and £ k , n L , 
and hm are given by equations (17). 

E. Temperature Distribution 

One of the assumptions made in developing the equations 
herein was that the material properties are not functions of 
temperature. This assumption was made only to conserve com- 
puter time. The assumption does, of course, restrict the 
maximum temperatures to around three hundred degrees Fahren- 
heit for aluminum. This range of temperature is, however, 
more than sufficient to demonstrate the validity of the 
theory before large deflection effects become significant. 

The program can handle either an anlytical temperature 
"surface" or an experimentally measured temperature distri- 
bution. The analytical temperature distribution is specified 
in the form of a polynomial in the independent space variables 
as shown in Fig. 2. The only requirement for the measured 
temperature is that the measurements be made at a sufficient 
number of points to accurately define the temperature distri- 
bution. 


The magnitude of the temperature distribution can be 
changed by input ing a series of AT's, In this case, since 
T re f is used as 1.0, the AT's are input as the actual value 
of the temperature desired (in degrees Fahrenheit or Centi- 
grade depending on the system of units used). 


Any experimentally measured temperature distribution may 
be input. The values of the temperature at the integration 
points are calculated by a two-dimensional, quadratic inter- 
polation subroutine. The temperatures • are input on a rectangu- 
lar grid. The points are evenly spaced in the £ and re- 
directions although the respective spacings need not be 
equal (i.e. the elements of the grid need not be square). 

A sample of the grid and an explanation of the defining 
parameters is shown in Fig. 3 • 


KC (I) 
LC( I) 

DTX 


- number of the first horizontal line at 
the Ith. vertical line (I = 1, 2, .., NTX) 

= the number of the last horizontal line 
at the vertical line (I =1, 2, .., 
NTX) 

= distance between vertical lines 


DTY 


distance between horizontal lines 
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(XT1, YT1) 


= coordinates of lower left hand point of 
the grid. 

NPTS = (not input) is calculated internally. 

This is the total number of grid points 

Fob the grid in Fig. 3, 


NTX 

= 

7 


KC( I) 

SB 

1, 2, 2, 3, 

3, 4, 4 

LK(I) 

= 

14, 14, 13, 

13, 12, 12, 

DTX 

= 

.142 


DTY 

-55 

.13 


( XT1 , YT1) 

= 

✓""N 

CO 

• 

\ 

<> 

Q 

• 

O 



The temperatures at the grid points are input from 
bottom-to-top for each vertical line starting from the 
left side. This sequence is shown by the circled numbers 
in Fig. 3. 


Interpolation will be attempted at any point within 
DTX and/or DTY of one of the grid points. This is a modi- 
fication of a program contained in Ref. 18, in which a 
complete description is given. 


The variable called TREF in the program is not actually 
used anywhere in the calculations. It is simply used as 
additional information to be output. Thus there are two 
ways of input ing the temperature distributions and incre- 
menting the AT values. 


The first method is to simply input the actual magni- 
tudes of the temperatures on the plate. In this case the 
values of AT will be of 0 (1). At AT=1° , then, the eigen- 
values calculated will give the frequencies of the plate for 
the input temperature distribution. 

If desired, the temperature distribution may be normal- 
ized with respect to the temperature, TREF, at some reference 
point on the plate. In this case as with the analytical 
distribution, the AT's are input as the actual value of the 
temperature desired at the reference point. 
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P. Thickness Distribution 


The thickness distribution h/h r , is symetric 
about the £-n plane and is described by two polynomials in 
6 and n . One gives the distribution on surface 1 and the 
other on surface 2. These two surfaces are separated by a _ 
line from the origin of the coordinate system at the angle 0. 
The value h (called TO in the program) is the thickness at 
the origin. 

G. Miscellaneous Comments 

The eigenvalue routine used here is a double precision 
version of the subroutine "NROOT" from the IBM Scientific 
Subroutine Package. (Note that this requires a double 
precision version of the subroutine "EIGEN” from the same 
source). The subroutine "DMINV" and "DGMPRD" (no listing 
given) are used directly from that source. 


Extensive use was made of the disk storage available 
in writing the program. This reduced the core storage re- 
quirements to around 250,000 bytes on the IBM 370-165 
computer used. Although the execution time for the program 
using all core storage would be about one-third of that 
using disk storage, the program would be limited to only 
thirty deflection and stress function terms and ten quadra- 
ture points. Also, core storage was a premium at the time 
of writing because of the large amount of business done by 
the Computer Center at The Ohio State University. 

It should also be noted that for the coordinate system 
used some plates will be symmetric about the x-axis. In 
these cases the even and odd terms in the assumed solution 
uncouple. Thus, the deflection function may be separated 
into one function containing only even terms in q and one 
containing only odd terms in q. Each of these functions 
can then be input separately to give all even modes or all 
odd modes respectively. The same comments also apply to the 
stress function. 


IV. RESULTS 


To compare to results in the literature, a conversion 
from the notation used in most other sources to that used 
herein is necessary. As long as the results are presented 
in non-dimensional form, only ratios of the material proper- 
ties are required. Thus let. 
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E 12 /E 22 * D xy /D 


2 D, /D 
k y 


All the data presented here will be converted using 
these relations to the notation previously described. 

All of the computations presented in this section were 
made using a 36 mixed term deflection function, 

5 5 

w U,n.) = g U.n) 2 e h. ,?V . 

i=o j=o 


Thus, the first thirty-six modes and frequencies were 
calculated. The runs took an average of three minutes 
(Central Processing Unit) time on the IBM 370-165- 

No effort was made to optimize the program, the purpose 
being to obtain consistantly good results for any planform 
shape with any boundary conditon. e.g,, acceptable results 
can be obtained for the torsion mode of a rectangular canti- 
lever plate with only three terms in the deflection function. 
However, this number of terms is completely inadequate for 
any other mode of the thermally stressed cantilever plate 
and is inadequate for any mode of any other of the many plates 
investigated. Thus, the large number of terms in both the 
stress function and the displacement function may be much 
greater than required for some of the problems solved. This 
point is immaterial when the choice boils down to either 
obtaining an accurate quantitative answer in which one can 
have confidence or some answer that may only be in the !, ball 
park." 

A. Comparison of Orthotropic Results Without Thermal Loading 

As was previously stated, the material published without 
thermal loading is voluminous. For the sake of brevity, only 
a few comparisons will be made . 

Tables 1 and 2 are comparisons of calculated frequencies 
from the literature with those calculated by this program. 

It can be seen that the method under discussion here gives 



excellent agreement with those frequencies. The expression 
for is given in the list of symbols and in Appendix A. 

Table 3 gives a comparison with some experimental 
frequencies for plywood plates. Note that the experimental 
values are higher than the calculated frequencies . Because 
the method used here gives solutions which converge from 
above the exact solution, these errors are attributed to 
restraints inherent in the experimental approximation of the 
simply supported boundary conditions. 

B. Comparison of Isotropic Results With Thermal Loading 

As was stated previously, no experimental data were 
found in the literature on the effect of thermal stresses 
on orthotropic plates. Thus, a brief quantitative compari- 
son is made with experimental data from Ref. 1 6 for the 
special case of the isotropic plate. The purpose is to show 
the agreement that was obtained for widely different cases. 
The isotropic plate elastic properties requires that, 

E 11 = E 22 = VU-v 4 ) 

E 12 = vE/U-v 2 ) 

G 12 = E/2( 1+v) 

A.s in Ref, 16, nominal values of the plate material proper- 
ties used were, 

E - ,10 7 psi 



a = 12.8x10 6 /°P 

Fig. 4 shows a comparison for the first five modes of 
a square cantilever plate. 16 It is interesting that the 
fourth and fifth mode frequency curves cross. A typical 
experimentally measured temperature distribution resulting 
from radiant lamp edge heating is shown in Fig. 5 • 

Fig. 6 presents an unsymetrical trapezoidal canti- 
lever plate that does not appear in Ref. 16. Only the 
first two modes were recorded for this plate. One of the 
temperature distributions measured on this plate is shown in 
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Pig. 7. 

Because of its boundary conditons and the choice of 
coordinate system, the plate shown in Pig. 8 is also unsymetri- 
cal. The stresses as well as the deflections are affected by 
the boundary conditions. 16 Here again, as in Pig. 4, two 
of the frequencies cross. The heating elements were centered 
over the diagonal from lower-left to upper-right giving a 
temperature distribution as shown in Pig. 9. 

The frequencies of a plate with a single point clamped 
is shown in Pig. 10. The agreement with the four modes 
measured is seen to be good. 

A plate with homogeneous pinned-free boundary condi- 
tions is shown in Pig. 11. Only the first two modes were 
recorded for this plate. The third calculated mode frequency 
is also shown. The temperature distribution shown in Pig. 

5 is also typical of that used to calculate the frequencies 
for the plates in Pigs. 10 and 11. 

C. Comparison of Orthotropic Results With Thermal Loading 

Pigs. 12, 13, 1^, 15 , 16 and 17 constitute the results 
of a very brief study that indicates the large effect that 
orthotropy can have in the presence of thermal gradients. 

For the sake of brevity, only one boundary condition, the 
cantilever plate, and only one assumed temperature distri- 
bution, T = AT j n j 3 , is presented. Also, only the two lowest 
modes are presented although as many of the higher modes as 
could possibly be desired are available in the computer 
print-out. It should be noted that orthotropy does not change 
the characteristic shape of the response curves shown. How- 
ever, because of the influence that the directional properties 
of the material can have on the stress field for a given 
temperature distribution, orthotropy can produce marked 
Increases In the loss of effective stiffness for a given 
heating rate. In each figure the isotropic response curves 
are given for comparison purposes. 

In Pigs. 12 and 13, it can be seen that doubling the 
thermal expansion coefficient, a 2 , in the chordwise direction 
has little effect on the frequency for this temperature dis- 
tribution because the chordwise stress component effect is 
small for this plate and remains small even when a 2 is 
doubled. However, when the longitudinal coefficient, ai , 
is doubled, the longitudinal stress component is essentially 
doubled and the frequency is seen to decrease at a markedly 
higher rate. Thermal buckling is reached at a temperature 
only one half as great as before. This means that, for a 
given heating rate, if an isotropic plate buckles within 
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thirty seconds, an orthotropic plate with this ratio of thermal 
coefficients would buckle in only fifteen seconds. 

In Fig. 14 it can be observed that doubling the chord- 
wise modulus of elasticity, E 22 , actually produces a decrease 
in the rate of frequency decay and an increase in .AT cr for the 
bending mode. That the various modes will not behave in the 
same way for a given material property change is shown in 
Fig. 15 where the doubled chordwise modulus produces the 
opposite effect on the torsion mode. Doubling the spanwise 
modulus, Epp, causes a higher rate of frequency decay with 
correspondingly lower buckling temperatures in both modes.. 

It should be noted, however, that, although large, doubling 
the longitudional modulus of elasticity does not have the 
extreme effect that Is caused by doubling the longitudinal 
conductivity coefficient. 

Fig. 16 shows that the decrease in the rate of decay 
of the bending mode, achieved by doubling E 22 in Fig. 14 , 
is more than offset by the increase caused by doubling 
in Fig. 12. Both Figs. 1 6 and 17 show appreciable destabili- 
zing effects when both the moduli and expansion coefficient 
are changed. It is recognized that the change in properties 
as used is large, but the effects are also large. Using 
the computer program presented herein, almost unlimited 
parametric studies may be made to determine trends or, more 
efficiently, calculations may be made to obtain answers for 
specific cases once a problem has been defined. 

D. Effect of Stress Distribution 

An interesting by-product of this program is an accurate 
calculation of the thermal stress distribution. The effect 
of the stress distribution function used in calculating the 
vibration frequency is shown in Fig. 15. The boundary condi- 
tions require a stress function containing both odd and even 
terms in n» Thus, the stress distribution calculated, using 
only even terms, does not give the correct stress distribution 
and hence, does not give the correct frequency. However, 
when the correct mixed terms were used, the correct stress 
distribution existing in the plate resulted and consequently , 
the calculated frequencies agreed more closely with the experi- 
mental values. The 36 mixed term results show that the 
analytically predicted stress distribution had effectively 
converged when 24 mixed terms were used. 

The three sets- of frequencies calculated were compared 
to the experimental response of the plate shown in Fig. 3. 

The results show that correct stresses are in fact necessary 
in order to obtain the correct frequencies . 
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V. CONCLUDING REMARKS 


The computer program presented herein is, in effect, a 
general solution to the problem of the linear vibration of 
thermally stressed trapezoidal plates. The theory has been 
verified experimentally for thermally stressed isotropic 
plates and has been found to agree favorably with analytical 
data found in the literature for orthotropic plates with no 
thermal loading. It appears that accurate results can be 
obtained by the methods herein described for almost any boundary 
conditions of practical importance. The solutions, based on 
linear theory, do not hold as the buckling region is approached 
because of the non-linear effects of large deflections. 

Experience in using the program shows that the number of 
terms required in the assumed solutions increases with the 
complexity of the geometry. However, consistently accurate 
results are obtained using a 30-36 term displacement function 
and a 24-30 term stress function. Accurate integration can 
be obtained for this number of terms using ten quadrature 
points in the ^-direction (twenty points in the n-direction) . 

A very extensive experimental program with orthotropic 
plates would be required to verify all the facets of the 
program as presented. However, the data comparisons shown 
combined with many cases of experimental isotropic plate 
data not presented herein gives the writers great confidence 
in the analytical results. 


The generality of this program should not be overlooked. 
Its extension to obtain the accurate solution of the flutter 
of thermally stressed plates and panels can be readily made-. 

A natural extension of this work would be to examine, without 
the assumption of mode identity, the large deflection effects 
observable as heating progresses, A recently developed 
method of solving large sets of non-linear equations shows 
great promise in the area of large deflections which has 
yet to be effectively investigated. 
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TABLE 1 


PF-CC-PP-CC 


AR 

E 11 

E 12 

E 22 

G 12 

h 

Ref. 

No. 

Calculated 

By 

Program 

Reference 

Values 

1.0 

1.0 

1.0 

1.0 

0.0 

28.93 

29.29 

8 

1.0 

2.0 

1.0 

1.0 

0 .0 

21.6 

21.82 

8 

1.0 

1.0 

3.0 

1.0 

0.0 

36.1 

36.5 

8 

1.0 

3.0 

3.0 

1.0 

0.0 

22.35 

22.56 

8 

1.0 

6 . o 

2.0 

1 . 0 

0.0 

1.6 , 12 

16.13 

8 

1.0 

1.0 

1.0 

2.0 ' 

0.0 

50.8 

51.5 

8 

1.0 

1.0 

3.0 

9.0 

0 . 0 

73.0 

74.0 

8 

2.0 

3.117 

0.12 

1.0 

B 

53.6 | 

53.7 

3 
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Comparison of Analytical and 
Experimental Response of a 
CC-FF-FF-FF Plate 
ct=-0.@ /3=0.0 7=0.0 
AR=5/3 , a -20", h=1/4" 

— measured 
© calculated 



Fig. 6 
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Typical Temperature 
Distribution used in 
Fig. 8 

T 

t80 



V 

Fig. 9 
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Comparison of Analytical and 
Experimental Response of a 

Plate Clamped at (0,0) 

AR=1.0 , a=18'\ h=3/16" 

— measured 
® calculated 



Plant orm 



Comparison of Analytical and 
Experimental Response of a 
PF-P F-PF-PF Plate 
AR=1.0 , a=1 8", h=3/16" 

— -measured 



AT ,°F 
Fig. n 


39 



Effect of a on Plate 
Vibration 
T=ATWl 3 

Constant Thickness 
-AR=5/3 

I s * mode 

En =E S2 



AT,°F 

Fig. 1.2 

4.0 




Effect of 'E' on Plate 
Vibration 
AR=5/3 

Constant Thickness 
T =ATI 7 ?I 3 

1 st mode E 12 =G 12 =0.375(10 7 )psi 

.^=02=12 8(1Cr®)/°F 

f Isotropic: 

" LEf]= E 22 =1.125 (10 7 )psi 



AT,°F 
Fig. 14 
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Effect of V E' on Plate 
Vibration 

2 nc * mode 

(See Fig 14 for notation) 
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Combined Effect of 'E' and a 
on Plate Vibration 


AR = 5/3 

Constant Thickness 
T =ATp7| 3 

I s * mode 

^ 2 = 0 , 2 = 0.375(10 7 )psi 
Isotropic: 

0 jE 11 = E 22 =1-125(1O 7 )psi 
0 U|= a 2 =12.8(10 6 ) l°F 



AT,°F 


Fig. 18 
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Combined Effect of 'E' and a 
on Plate Vibration 

2 nd mode 

( See Fig. 16 for notation ) 
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Effect of Stress Function 
on Plate Vibration 

measured 

calculated using: 
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Appendix A 


Matrix Elements and Parameters 


3 


E, 


B ij ,k£ ~ //< h r > (c Ws? + V ( a k«. 5 nn 


+ <!) 2 ^ t(a ±J ) K <«„) +<.„> (^,M 


+ 4( l> 2 <“ij>5n Kl’s,’ dWn 


M ij ,kz //{P nn ( ‘ a ij h (o Wg + ^ a ij ^ n (a kj2 n 


(c Wn + (a kJlV } 


T ij,kil ff h r (a ij 5 ^ a k£ ^ 


A pq,rs = ff “FT {( |^ ^rvi)nn(Y«J nr, + ^ < Y„ n ) f ? < Y^* ) i 


r pq nn rs Tin a-^ ’pq rs g : g 


a 2 

+ a- £ ! " J ^^pq^E^rPnn + ' v ::*pnn '"'rd'iC' 
+ *11 ‘Vl^n^rsW aWn 


a. 


r rs ■ "d|) (y rs ) m + ^ (Y rs ) 55 ] T <«,„) d 5 dn 
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k 2 = ( a l AT/a 1]L ) 

{C}= -A— {€} 

a\ 

a ij “ sU,n) 

Ypq = f(€,n)S p n q 
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Appendix B 


Logic Plow Diagram 


C START J 

I 


' READ NUMBER 
OP DIFFERENT 
BOUNDARY CONDITIONS 


DO THRU a 
FOR EACH SET OF 
BOUNDARY CONDITIONS 


READ PARAMETERS FOR 
BOUNDARY CONDITION AND 
EQUILIBRATING FUNCTIONS 


I 


READ a, 0, y, bp t Q ; 
INTEGRATION DATA, AND 
THICKNESS DATA 


READ NUMBER OF DIFFERENT 
SETS OF ASPECT RATIOS AND/OR 
MATERIAL PROPERTIES 
AND THEIR VALUES 


READ NUMBER OF DEFLECTION 
FUNCTION TERMS, AND £ AND n 
EXPONENTS 


T 
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DO THRU Y FOR 
EACH AT VALUE 
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Appendix C 


Program Listing 
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